Rasch analysis of outcomes

1 Setup

Code

library(easyRasch) # devtools::install_github("pgmj/easyRasch")
#library(RISEkbmRasch)
library(grateful)
library(ggrepel)
library(car)
library(kableExtra)
library(tidyverse)
library(eRm)
library(iarm)
library(mirt)
library(psych)
library(psychotree)
library(matrixStats)
library(reshape)
library(knitr)
library(patchwork)
library(formattable) 
library(glue)


### optional libraries
#library(TAM)
#library(skimr)
#library(janitor)
# 1727699
# nohup quarto render PDSSSR.qmd --to html &> render_pdsssr.out & 

### some commands exist in multiple packages, here we define preferred ones that are frequently used
select <- dplyr::select
count <- dplyr::count
recode <- car::recode
rename <- dplyr::rename

path_prim <- '/srv/projects/nils/study4rasch'
path2 <- '/Volumes/projects/k8_CPF_Kaldo/Data/Lärande Maskiner/Nils/Study4data'
running_machine = Sys.info()[['sysname']]
path = ifelse(running_machine!='Linux',path2,path_prim)
set_mywd_path = ifelse(running_machine!='Linux','./study4rasch/rasch_measures','./rasch_measures')
n_cores =ifelse(running_machine!='Linux',1,8)

tryCatch({setwd(set_mywd_path)}, error = function(set_mywd_path){print('already in right directory')})

[1] "already in right directory"

Code

source('settings.R') # containing item labels 
source('mod_easyRasch_func.R') # modified functions + added 

2 Importing data

Code

### import data
df_raw <- read.csv(file=file.path(path,'data','raw_data','outcome.csv')) #

### Clean data 
source("preprocess.R")
df <- preprocess_data(df_raw)

# items to use 
items_to_use <- 'PDSS.SR' # options are #MADRS.S, LSAS.SR, PDSS.SR
itemlabels <- itemlabels %>% 
  as_tibble() %>%
  filter(str_detect(itemnr, items_to_use))


# DIF+aux variables to use 
# DIF not used 'marital_status','children', 'working'
dif_variables <- c('Patient','Treatment','sex','age','TreatmentAccessStart','education')

### Make a backup of the dataframe, in case you need to revert changes at some point
df.all <- df

print(itemlabels)

# A tibble: 7 × 2
  itemnr    item                
  <chr>     <chr>               
1 PDSS.SR_1 Number              
2 PDSS.SR_2 Distress            
3 PDSS.SR_3 Worry/Fear          
4 PDSS.SR_4 Avoidance places    
5 PDSS.SR_5 Avoidance activity  
6 PDSS.SR_6 Impairment work/home
7 PDSS.SR_7 Impairment social   

3 Checking missing

Code

itemStart <- items_to_use
out = RImissing_mod(df,itemStart,facet="Time")
print(out)

As we see here PDSS.SR is not used for all conditions over the course of their treatments. However during screening they are. We start by showing all timepoints, to draw targeting plots.

3.1 Remove missing

Code

##### Before filtering out participants, you should check the missing data structure using RImissing() and RImissingP()

# If you want to include participants with missing data, input the minimum number of items responses that a participant should have to be included in the analysis:
min.responses <- 4 # In our case if they have missing on one item all items for that time are missing. 

df_save <- df
# Select the variables we will work with, and filter out respondents with missing data
df <- df %>% 
  select(all_of(c(dif_variables,"Time",itemlabels$itemnr))) %>%  # v
  filter(rowSums(is.na(select(.,all_of(itemlabels$itemnr)))) < min.responses) 

4 Create DIF variables

Code

#---- Create DIF variables----
  
# DIF variables into vectors, recoded as factors since DIF functions need this
# these could also be stored in its own dataframe (not a tibble) instead of as vectors
# Named vector for the new types

type_transform <- c(Treatment = "factor", sex = "factor", age = "numeric",
                    TreatmentAccessStart ="integer",education="factor",Time='factor')

# Transform columns based on the named vector
dif <- df %>%
  mutate(across(names(type_transform), ~ switch(type_transform[which(names(type_transform) == cur_column())], 
                                               "integer" = as.integer(.), 
                                               "character" = as.character(.),
                                               "factor" = as.factor(.), 
                                               "numeric" = as.numeric(.)))) %>%
                                               as.data.frame(.) %>%
                                               select(!all_of(itemlabels$itemnr))


# then remove them from dataframe, since we need a dataframe with only item data for the Rasch analyses
df <- df %>% select(all_of(c("Time",itemlabels$itemnr))) # add time here ? 
dfnotime <- df %>% select(!Time)
source("RISE_theme.R")

5 Demographics descriptives

Code

dif_spec <- dif %>% filter(Time=='SCREEN') %>% select(!all_of(c("Time","Patient")))
summary(dif_spec)
# RIdemographics(dif_spec,'Treatment') <- function crashes TODO make own 

 Treatment  sex           age        TreatmentAccessStart         education   
 MDD:2853   F:3834   Min.   :18.00   Min.   :2008         Primary      : 484  
 PD :1592   M:2366   1st Qu.:27.00   1st Qu.:2012         Secondary    :3024  
 SAD:1755            Median :33.00   Median :2014         Postsecondary:2692  
                     Mean   :35.24   Mean   :2014                             
                     3rd Qu.:42.00   3rd Qu.:2017                             
                     Max.   :84.00   Max.   :2019                             

6 Overall number of responses

Code

# Collapsed 
RIallresp(dfnotime)
# Seperate 
RIallresp_over_times <- df %>% # 
  split(.$Time) %>% # split the data 
  map(~ RIallresp(.x %>% dplyr::select(!Time))) #+ labs(title = .x$Time)) # create separate for each time

# make nice later 
RI_allresp_kable_grid = combine_kables_grid(RIallresp_over_times,cols=3)
RI_allresp_kable_grid

Response category	Number of responses	Percent
0	48333	31.8
1	54744	36.1
2	34135	22.5
3	12359	8.1
4	2182	1.4

x

SCREEN.Response category SCREEN.Number of responses SCREEN.Percent PRE.Response category PRE.Number of responses PRE.Percent WEEK01.Response category WEEK01.Number of responses WEEK01.Percent 0 16055 37.0 0 1661 14.1 0 1915 19.8 1 10996 25.3 1 3725 31.7 1 3879 40.1 2 10206 23.5 2 4114 35.0 2 2868 29.6 3 5057 11.7 3 1921 16.4 3 904 9.3 4 1086 2.5 4 325 2.8 4 115 1.2 WEEK02.Response category WEEK02.Number of responses WEEK02.Percent WEEK03.Response category WEEK03.Number of responses WEEK03.Percent WEEK04.Response category WEEK04.Number of responses WEEK04.Percent 0 2114 21.6 0 2413 24.7 0 2649 28.0 1 4268 43.6 1 4262 43.6 1 4041 42.8 2 2545 26.0 2 2355 24.1 2 2093 22.1 3 768 7.8 3 665 6.8 3 573 6.1 4 98 1.0 4 70 0.7 4 94 1.0 WEEK05.Response category WEEK05.Number of responses WEEK05.Percent WEEK06.Response category WEEK06.Number of responses WEEK06.Percent WEEK07.Response category WEEK07.Number of responses WEEK07.Percent 0 2679 29.6 0 2903 33.3 0 2841 34.1 1 3911 43.2 1 3656 41.9 1 3482 41.8 2 1887 20.8 2 1669 19.1 2 1540 18.5 3 487 5.4 3 425 4.9 3 397 4.8 4 87 1.0 4 76 0.9 4 70 0.8 WEEK08.Response category WEEK08.Number of responses WEEK08.Percent WEEK09.Response category WEEK09.Number of responses WEEK09.Percent WEEK10.Response category WEEK10.Number of responses WEEK10.Percent 0 3020 37.8 0 2962 38.8 0 3027 42.3 1 3273 40.9 1 3082 40.4 1 2842 39.7 2 1348 16.9 2 1230 16.1 2 1042 14.6 3 313 3.9 3 325 4.3 3 220 3.1 4 40 0.5 4 31 0.4 4 30 0.4 POST.Response category POST.Number of responses POST.Percent " " " " 0 4094 45.4 1 3327 36.9 2 1238 13.7 3 304 3.4 4 60 0.7

7 Descriptives of raw data

Response distribution for all items are summarized below.

7.1 Floor/ceiling effects

Code

RIrawdist(dfnotime)

Code

#df_test = df %>% filter(Time=='SCREEN') %>% select(!Time)
# Across time 

# seperated time - might need to fix this remake in ggplot?
png(filename=paste0("./plots/",itemStart,"_rawdist_over_time.png"),
width=800,height=2400)
par(mfrow=c(7,2),mar = c(4, 4, 4, 2))  # c(5 <- bottom lines, 4, 4, 2)
RIrawdist_over_times <- df %>% # 
  split(.$Time) %>% # split the data 
  imap(~ {
    RIrawdist(.x %>% dplyr::select(!Time))
    mtext(paste("Time:", .y), side = 1, line = 3, adj = 0, cex = 0.8)
    })

dev.off()

png 
  2 

Code

print('Saved plot to /plots')

[1] "Saved plot to /plots"

7.2 Guttman structure

Code

RIheatmap(dfnotime)

Code

ggplots_with_spec_func_over_time(df,func_call = RIheatmap,title_all='Guttman structure over time',blankx=TRUE)

7.3 Descriptives - item level

Code

df_test = df %>% filter(Time=='SCREEN') %>% select(!Time)
RIlistItemsMargin(df, fontsize = 12)

itemnr	item
PDSS.SR_1	Number
PDSS.SR_2	Distress
PDSS.SR_3	Worry/Fear
PDSS.SR_4	Avoidance places
PDSS.SR_5	Avoidance activity
PDSS.SR_6	Impairment work/home
PDSS.SR_7	Impairment social

Tile plot

Code

# over all times 
RItileplot(dfnotime)

Code

# by time 
ggplots_with_spec_func_over_time(df,func_call = RItileplot,title_all='',blankx=FALSE)

Stacked bars

Code

RIbarstack(dfnotime)

Code

# by time 
ggplots_with_spec_func_over_time(df,func_call = RIbarstack,title_all='',blankx=FALSE)

Barplots

Code

RIbarplot(dfnotime)

8 Rasch analysis 1

Code

df_test = df %>% filter(Time=='SCREEN') %>% select(!Time)
RIlistItemsMargin(df, fontsize = 12)

itemnr	item
PDSS.SR_1	Number
PDSS.SR_2	Distress
PDSS.SR_3	Worry/Fear
PDSS.SR_4	Avoidance places
PDSS.SR_5	Avoidance activity
PDSS.SR_6	Impairment work/home
PDSS.SR_7	Impairment social

The eRm package, which uses Conditional Maximum Likelihood (CML) estimation, will be used primarily. For this analysis, the Partial Credit Model will be used.

8.1 Timepoint decision

Due to the amount of data across timepoints targeting will first be inspected, and the best fitting timepoint will be chosen based on targeting, which in turn will inform the rest of the analyses.

8.1.1 Targeting

Code

# increase fig-height above as needed, if you have many items og. value was 5. 
#RItargeting(dfnotime)

require(patchwork)
targeting_time_plots <- ggplots_with_spec_func_over_time(
  df,func_call = RItargeting,title_all='',
  blankx=FALSE,colsplot = 3,rowsplot=5,return_list_plots=TRUE)
good_layout_targeting_plots <- wrap_plots(targeting_time_plots, ncol = 3) + plot_layout(guides= "collect")
#print(good_layout_targeting_plots)
ggsave(
  filename = paste0("./plots/",items_to_use,"_targeting_over_time.png"),
  plot = good_layout_targeting_plots,
  width = 20,  # Increase width (in inches)
  height = 45, # Increase height (in inches)
  dpi = 300    # High resolution
)
print(good_layout_targeting_plots)

Pre has very good targeting looks promising. We also note that this is only for patient in treatment for panic disorder (no other treatment-conditions are mixed in PRE). We doublecheck how much data we have for each threshold. (Recomm. 20+)

8.1.2 Number of responses per threshhold

Code

dataframen4 <- df %>% filter(Time=="PRE") %>% select(!Time)
sapply(dataframen4, function(column) table(column))

  PDSS.SR_1 PDSS.SR_2 PDSS.SR_3 PDSS.SR_4 PDSS.SR_5 PDSS.SR_6 PDSS.SR_7
0       129       157        52       161       521       368       273
1       673       389       409       504       526       621       603
2       679       807       628       660       430       426       484
3       148       296       496       293       178       227       283
4        49        29        93        60        23        36        35

Based on the above we move forward with the PRE timepoint.

8.2 Analysis continuation PRE

Code

df_pre <- df %>% filter(Time=='PRE') %>% select(!Time) #%>% sample_n(800)

Conditional item fit

Code

RIitemfit(df_pre, cutoff = "Smith98")

Item	InfitMSQ	Location	1 ± 2 / √n
PDSS.SR_1	0.996	-0.10	[0.951, 1.049]
PDSS.SR_2	0.875	-0.07	[0.951, 1.049]
PDSS.SR_3	0.973	-0.96	[0.951, 1.049]
PDSS.SR_4	1.045	-0.24	[0.951, 1.049]
PDSS.SR_5	1.29	0.75	[0.951, 1.049]
PDSS.SR_6	0.956	0.40	[0.951, 1.049]
PDSS.SR_7	0.915	0.22	[0.951, 1.049]
Note:
MSQ values based on conditional estimation (n = 1678 complete cases).

Code

simfit1 <- RIgetfit(df_pre, iterations = 300, cpu = n_cores)  

RIitemfit(df_pre, simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
PDSS.SR_1	0.996	[0.929, 1.081]	0.981	[0.927, 1.075]	no misfit	no misfit	-0.10
PDSS.SR_2	0.875	[0.921, 1.085]	0.874	[0.918, 1.091]	0.046	0.044	-0.07
PDSS.SR_3	0.973	[0.946, 1.06]	0.976	[0.946, 1.065]	no misfit	no misfit	-0.96
PDSS.SR_4	1.045	[0.919, 1.084]	1.057	[0.914, 1.088]	no misfit	no misfit	-0.24
PDSS.SR_5	1.29	[0.899, 1.093]	1.342	[0.875, 1.131]	0.197	0.211	0.75
PDSS.SR_6	0.956	[0.926, 1.074]	0.936	[0.905, 1.087]	no misfit	no misfit	0.40
PDSS.SR_7	0.915	[0.932, 1.085]	0.916	[0.921, 1.08]	0.017	0.005	0.22
Note:
MSQ values based on conditional calculations (n = 1678 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RIgetfitPlot(simfit1, df_pre)

CICC

Code

ICCplot(as.data.frame(df_pre), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[1:4])

[1] "Please press Zoom on the Plots window to see the plot"

Code

### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)

Code

ICCplot(as.data.frame(df_pre), 
        itemnumber = 5:7, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[5:7])

[1] "Please press Zoom on the Plots window to see the plot"

Item-restscore

Code

RIrestscore(df_pre)

Item	Observed value	Model expected value	Absolute difference	Adjusted p-value (BH)	Statistical significance level	Location	Relative location
PDSS.SR_1	0.63	0.58	0.05	0.011	*	-0.10	0.62
PDSS.SR_2	0.66	0.59	0.07	0.000	***	-0.07	0.64
PDSS.SR_3	0.60	0.59	0.01	0.311		-0.96	-0.25
PDSS.SR_4	0.58	0.59	0.01	0.715		-0.24	0.47
PDSS.SR_5	0.47	0.60	0.13	0.000	***	0.75	1.46
PDSS.SR_6	0.62	0.59	0.03	0.116		0.40	1.11
PDSS.SR_7	0.64	0.59	0.05	0.011	*	0.22	0.93

Item-restscore-bootstrapped

Code

RIbootRestscore(df_pre,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600

Item	Item-restscore result	% of iterations	Conditional MSQ infit	Relative average item location
PDSS.SR_1	no misfit	62.8	1.00	0.62
PDSS.SR_1	overfit	37.2	1.00	0.62
PDSS.SR_2	no misfit	16.6	0.87	0.64
PDSS.SR_2	overfit	83.4	0.87	0.64
PDSS.SR_3	no misfit	92.6	0.97	-0.25
PDSS.SR_3	overfit	7.4	0.97	-0.25
PDSS.SR_4	no misfit	99.4	1.05	0.47
PDSS.SR_5	underfit	99.8	1.29	1.46
PDSS.SR_6	no misfit	85.6	0.96	1.11
PDSS.SR_6	overfit	14.4	0.96	1.11
PDSS.SR_7	no misfit	63.2	0.91	0.93
PDSS.SR_7	overfit	36.8	0.91	0.93
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 1678.

PCA

Code

RIpcmPCA(df_pre)

PCA of Rasch model residuals

Eigenvalues	Proportion of variance
1.90	26.1%
1.28	21.3%
1.21	17.1%
1.07	14.7%
0.89	12.2%

Residual correlations

Code

simcor1 <- RIgetResidCor(df_pre, iterations = 500, cpu = n_cores)
RIresidcorr(df_pre, cutoff = simcor1$p99)

	PDSS.SR_1	PDSS.SR_2	PDSS.SR_3	PDSS.SR_4	PDSS.SR_5	PDSS.SR_6	PDSS.SR_7
PDSS.SR_1
PDSS.SR_2	0.37
PDSS.SR_3	-0.08	-0.05
PDSS.SR_4	-0.32	-0.28	-0.08
PDSS.SR_5	-0.24	-0.23	-0.22	-0.15
PDSS.SR_6	-0.15	-0.17	-0.2	-0.17	-0.14
PDSS.SR_7	-0.26	-0.26	-0.16	0.09	-0.2	-0.03
Note:
Relative cut-off value is -0.046, which is 0.094 above the average correlation (-0.14). Correlations above the cut-off are highlighted in red text.

Partial gamma LD

Code

RIpartgamLD(df_pre)

Item 1	Item 2	Partial gamma	SE	Lower CI	Upper CI	Adjusted p-value (BH)
PDSS.SR_2	PDSS.SR_1	0.737	0.026	0.686	0.787	0.000
PDSS.SR_1	PDSS.SR_2	0.724	0.027	0.670	0.778	0.000
PDSS.SR_7	PDSS.SR_4	0.34	0.038	0.266	0.415	0.000
PDSS.SR_4	PDSS.SR_7	0.305	0.038	0.230	0.380	0.000
PDSS.SR_7	PDSS.SR_6	0.212	0.039	0.135	0.288	0.000
PDSS.SR_6	PDSS.SR_7	0.187	0.040	0.108	0.265	0.000
PDSS.SR_2	PDSS.SR_3	0.164	0.043	0.080	0.248	0.006
PDSS.SR_3	PDSS.SR_4	0.141	0.041	0.062	0.221	0.021

1st contrast loadings

Code

RIloadLoc(df_pre)

Response categories

Code

mirt(df_pre, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code

# for fewer items or a more magnified figure, use:
#RIitemCats(df)

Targeting

Code

# increase fig-height above as needed, if you have many items
RItargeting(df_pre)

Item hierarchy

Code

RIitemHierarchy(df_pre)

Score groups

Code

iarm::score_groups(as.data.frame(df_pre)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)

  score_group   n
1           1 894
2           2 784

Code

dif_plots <- df_pre %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code

clr_tests(df_pre, model = "PCM")

Conditional Likelihood Ratio Tests:

        clr df pvalue  sig 
overall 115 27 7.4e-13  ***

Code

item_obsexp(PCM(df_pre))

Score group 1: 
          mean obs mean exp std.res sig
PDSS.SR_1  1.1466   1.1648  -0.9382    
PDSS.SR_2  1.3078   1.3388  -1.4004    
PDSS.SR_3  1.5896   1.5887   0.0406    
PDSS.SR_4  1.2255   1.2332  -0.3487    
PDSS.SR_5  0.6742   0.5928   3.8498 ++ 
PDSS.SR_6  0.7835   0.7779   0.2626    
PDSS.SR_7  0.9267   0.9450  -0.8361    

Score group 2: 
          mean obs mean exp std.res sig
PDSS.SR_1  2.107    2.087    0.882     
PDSS.SR_2  2.354    2.317    1.701     
PDSS.SR_3  2.696    2.685    0.466     
PDSS.SR_4  2.365    2.357    0.360     
PDSS.SR_5  1.801    1.893   -3.443  -- 
PDSS.SR_6  2.042    2.048   -0.238     
PDSS.SR_7  2.215    2.193    0.841     

Code

grouping_based_on_score = score_groups(as.data.frame(df_pre))
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = grouping_based_on_score) 

       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 PDSS.SR_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 PDSS.SR_2 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 PDSS.SR_3 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 PDSS.SR_4 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
5 PDSS.SR_5 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
6 PDSS.SR_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
7 PDSS.SR_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA

DIF-sex

Code

dif.sex <- dif %>% filter(Time=='PRE') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_pre, dif.sex)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.sex) 

       Item     Var   gamma     se pvalue padj.BH sig   lower   upper
1 PDSS.SR_1 dif.sex  0.0470 0.0522 0.3676  1.0000     -0.0553  0.1494
2 PDSS.SR_2 dif.sex  0.0564 0.0519 0.2777  1.0000     -0.0454  0.1582
3 PDSS.SR_3 dif.sex -0.0701 0.0489 0.1519  1.0000     -0.1661  0.0258
4 PDSS.SR_4 dif.sex -0.1046 0.0487 0.0316  0.2215     -0.2001 -0.0092
5 PDSS.SR_5 dif.sex  0.0693 0.0456 0.1280  0.8959     -0.0199  0.1586
6 PDSS.SR_6 dif.sex  0.0318 0.0483 0.5105  1.0000     -0.0629  0.1265
7 PDSS.SR_7 dif.sex -0.0061 0.0490 0.9015  1.0000     -0.1020  0.0899

DIF-age

Code

dif.age <- dif %>% filter(Time=='PRE') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_pre, dif.age)

Item	2	3	Mean location	StDev	MaxDiff
PDSS.SR_1	-0.138	-0.001	-0.069	0.096	0.136
PDSS.SR_2	-0.108	0.070	-0.019	0.125	0.177
PDSS.SR_3	-1.122	-0.649	-0.885	0.335	0.474
PDSS.SR_4	-0.164	-0.445	-0.304	0.199	0.282
PDSS.SR_5	0.792	0.659	0.726	0.094	0.133
PDSS.SR_6	0.472	0.243	0.357	0.162	0.229
PDSS.SR_7	0.267	0.123	0.195	0.102	0.144

Code

partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.age) 

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 PDSS.SR_1 dif.age -0.0639 0.0311 0.0402  0.2813      -0.1248 -0.0029
2 PDSS.SR_2 dif.age -0.0420 0.0305 0.1693  1.0000      -0.1018  0.0179
3 PDSS.SR_3 dif.age -0.1140 0.0278 0.0000  0.0003  *** -0.1686 -0.0595
4 PDSS.SR_4 dif.age  0.0919 0.0291 0.0016  0.0109    *  0.0350  0.1489
5 PDSS.SR_5 dif.age  0.0096 0.0277 0.7281  1.0000      -0.0447  0.0640
6 PDSS.SR_6 dif.age  0.0401 0.0277 0.1476  1.0000      -0.0142  0.0943
7 PDSS.SR_7 dif.age  0.0438 0.0287 0.1260  0.8817      -0.0123  0.1000

DIF-education

Code

dif.edu <- dif %>% filter(Time=='PRE') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_pre, dif.edu)

Item	2	3	Mean location	StDev	MaxDiff
PDSS.SR_1	-0.078	-0.192	-0.135	0.081	0.114
PDSS.SR_2	-0.116	0.200	0.042	0.223	0.316
PDSS.SR_3	-0.908	-1.089	-0.999	0.128	0.181
PDSS.SR_4	-0.236	-0.304	-0.270	0.048	0.068
PDSS.SR_5	0.769	0.768	0.768	0.001	0.001
PDSS.SR_6	0.396	0.383	0.389	0.009	0.013
PDSS.SR_7	0.172	0.234	0.203	0.043	0.061

Code

partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.edu) 

       Item     Var   gamma     se pvalue padj.BH sig   lower   upper
1 PDSS.SR_1 dif.edu  0.0646 0.0470 0.1688  1.0000     -0.0274  0.1567
2 PDSS.SR_2 dif.edu  0.0282 0.0465 0.5443  1.0000     -0.0630  0.1195
3 PDSS.SR_3 dif.edu  0.0292 0.0435 0.5025  1.0000     -0.0561  0.1144
4 PDSS.SR_4 dif.edu -0.0733 0.0432 0.0900  0.6297     -0.1580  0.0114
5 PDSS.SR_5 dif.edu  0.1041 0.0410 0.0112  0.0783   .  0.0237  0.1844
6 PDSS.SR_6 dif.edu -0.0347 0.0433 0.4225  1.0000     -0.1196  0.0502
7 PDSS.SR_7 dif.edu -0.1511 0.0431 0.0005  0.0032  ** -0.2357 -0.0666

DIF-Year

Code

dif.yr <- dif %>% filter(Time=='PRE') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_pre, dif.yr)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.yr) 

       Item    Var   gamma     se pvalue padj.BH sig   lower  upper
1 PDSS.SR_1 dif.yr  0.0317 0.0329 0.3354  1.0000     -0.0328 0.0963
2 PDSS.SR_2 dif.yr -0.0079 0.0334 0.8120  1.0000     -0.0733 0.0574
3 PDSS.SR_3 dif.yr  0.0320 0.0296 0.2796  1.0000     -0.0260 0.0899
4 PDSS.SR_4 dif.yr -0.0377 0.0308 0.2206  1.0000     -0.0979 0.0226
5 PDSS.SR_5 dif.yr -0.0440 0.0277 0.1117  0.7819     -0.0983 0.0102
6 PDSS.SR_6 dif.yr  0.0475 0.0306 0.1206  0.8442     -0.0125 0.1075
7 PDSS.SR_7 dif.yr -0.0204 0.0301 0.4984  1.0000     -0.0795 0.0387

8.3 Analysis part 1 decision

From conditional item fit item 2 and 5 are flagged. However Item 2 have large residual correlations with 1 (possibly also seen at the contrast loadings), however item 5 is most different from the score groups. Removing item 5 due to the underfit.

9 Rasch analysis 2

Based on the above we remove item 5

Code

df_pre <- df %>% filter(Time=='PRE') %>% select(!c(Time,PDSS.SR_5)) #%>% sample_n(800)
RIlistItemsMargin(df_pre, fontsize = 12)

itemnr	item
PDSS.SR_1	Number
PDSS.SR_2	Distress
PDSS.SR_3	Worry/Fear
PDSS.SR_4	Avoidance places
PDSS.SR_6	Impairment work/home
PDSS.SR_7	Impairment social

Conditional item fit

Code

RIitemfit(df_pre, cutoff = "Smith98")

Item	InfitMSQ	Location	1 ± 2 / √n
PDSS.SR_1	1.015	0.04	[0.951, 1.049]
PDSS.SR_2	0.897	0.04	[0.951, 1.049]
PDSS.SR_3	1.013	-0.88	[0.951, 1.049]
PDSS.SR_4	1.12	-0.12	[0.951, 1.049]
PDSS.SR_6	1.043	0.56	[0.951, 1.049]
PDSS.SR_7	0.961	0.36	[0.951, 1.049]
Note:
MSQ values based on conditional estimation (n = 1678 complete cases).

Code

simfit1 <- RIgetfit(df_pre, iterations = 300, cpu = n_cores)  

RIitemfit(df_pre, simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
PDSS.SR_1	1.015	[0.93, 1.074]	0.992	[0.929, 1.07]	no misfit	no misfit	0.04
PDSS.SR_2	0.897	[0.922, 1.069]	0.896	[0.92, 1.071]	0.025	0.024	0.04
PDSS.SR_3	1.013	[0.931, 1.074]	1.012	[0.93, 1.079]	no misfit	no misfit	-0.88
PDSS.SR_4	1.12	[0.93, 1.071]	1.132	[0.93, 1.071]	0.049	0.061	-0.12
PDSS.SR_6	1.043	[0.925, 1.096]	1.017	[0.918, 1.107]	no misfit	no misfit	0.56
PDSS.SR_7	0.961	[0.931, 1.079]	0.968	[0.919, 1.089]	no misfit	no misfit	0.36
Note:
MSQ values based on conditional calculations (n = 1678 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RIgetfitPlot(simfit1, df_pre)

CICC

Code

ICCplot(as.data.frame(df_pre), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[1:4])

[1] "Please press Zoom on the Plots window to see the plot"

Code

### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)

Code

ICCplot(as.data.frame(df_pre), 
        itemnumber = 5:6, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[5:6])

[1] "Please press Zoom on the Plots window to see the plot"

Item-restscore

Code

RIrestscore(df_pre)

Item	Observed value	Model expected value	Absolute difference	Adjusted p-value (BH)	Statistical significance level	Location	Relative location
PDSS.SR_1	0.64	0.61	0.03	0.189		0.04	0.65
PDSS.SR_2	0.67	0.62	0.05	0.003	**	0.04	0.66
PDSS.SR_3	0.61	0.61	0.00	0.868		-0.88	-0.26
PDSS.SR_4	0.58	0.62	0.04	0.117		-0.12	0.49
PDSS.SR_6	0.61	0.62	0.01	0.697		0.56	1.17
PDSS.SR_7	0.64	0.62	0.02	0.282		0.36	0.97

Item-restscore-bootstrapped

Code

RIbootRestscore(df_pre,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600

Item	Item-restscore result	% of iterations	Conditional MSQ infit	Relative average item location
PDSS.SR_1	no misfit	94.8	1.02	0.65
PDSS.SR_1	overfit	5.2	1.02	0.65
PDSS.SR_2	no misfit	52.8	0.90	0.66
PDSS.SR_2	overfit	47.2	0.90	0.66
PDSS.SR_3	no misfit	99.4	1.01	-0.26
PDSS.SR_4	no misfit	91.0	1.12	0.49
PDSS.SR_4	underfit	9.0	1.12	0.49
PDSS.SR_6	no misfit	99.8	1.04	1.17
PDSS.SR_7	no misfit	96.4	0.96	0.97
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 1678.

PCA

Code

RIpcmPCA(df_pre)

PCA of Rasch model residuals

Eigenvalues	Proportion of variance
1.96	32.4%
1.30	22.1%
1.11	18.8%
0.93	15.8%
0.68	10.7%

Residual correlations

Code

simcor1 <- RIgetResidCor(df_pre, iterations = 500, cpu = n_cores)
RIresidcorr(df_pre, cutoff = simcor1$p99)

	PDSS.SR_1	PDSS.SR_2	PDSS.SR_3	PDSS.SR_4	PDSS.SR_6	PDSS.SR_7
PDSS.SR_1
PDSS.SR_2	0.33
PDSS.SR_3	-0.13	-0.11
PDSS.SR_4	-0.38	-0.34	-0.12
PDSS.SR_6	-0.2	-0.22	-0.25	-0.2
PDSS.SR_7	-0.33	-0.32	-0.22	0.05	-0.07
Note:
Relative cut-off value is -0.085, which is 0.083 above the average correlation (-0.167). Correlations above the cut-off are highlighted in red text.

Partial gamma LD

Code

RIpartgamLD(df_pre)

Item 1	Item 2	Partial gamma	SE	Lower CI	Upper CI	Adjusted p-value (BH)
PDSS.SR_2	PDSS.SR_1	0.737	0.026	0.686	0.788	0
PDSS.SR_1	PDSS.SR_2	0.707	0.028	0.652	0.763	0
PDSS.SR_7	PDSS.SR_4	0.352	0.037	0.279	0.426	0
PDSS.SR_4	PDSS.SR_7	0.303	0.038	0.228	0.379	0
PDSS.SR_7	PDSS.SR_6	0.235	0.039	0.158	0.312	0
PDSS.SR_6	PDSS.SR_7	0.179	0.041	0.100	0.259	0

1st contrast loadings

Code

RIloadLoc(df_pre)

Response categories

Code

mirt(df_pre, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code

# for fewer items or a more magnified figure, use:
#RIitemCats(df)

Targeting

Code

# increase fig-height above as needed, if you have many items
RItargeting(df_pre)

Item hierarchy

Code

RIitemHierarchy(df_pre)

Score groups

Code

iarm::score_groups(as.data.frame(df_pre)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)

  score_group   n
1           1 913
2           2 765

Code

dif_plots <- df_pre %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code

clr_tests(df_pre, model = "PCM")

Conditional Likelihood Ratio Tests:

        clr df pvalue  sig 
overall 60  23 3.8e-05  ***

Code

item_obsexp(PCM(df_pre))

Score group 1: 
          mean obs mean exp std.res sig
PDSS.SR_1  1.151    1.154   -0.171     
PDSS.SR_2  1.319    1.333   -0.658     
PDSS.SR_3  1.600    1.581    0.944     
PDSS.SR_4  1.236    1.224    0.597     
PDSS.SR_6  0.779    0.766    0.665     
PDSS.SR_7  0.914    0.934   -0.952     

Score group 2: 
          mean obs mean exp std.res sig
PDSS.SR_1  2.123    2.119    0.164     
PDSS.SR_2  2.364    2.347    0.828     
PDSS.SR_3  2.709    2.723   -0.639     
PDSS.SR_4  2.379    2.394   -0.629     
PDSS.SR_6  2.076    2.092   -0.618     
PDSS.SR_7  2.260    2.236    0.961     

Code

grouping_based_on_score = score_groups(as.data.frame(df_pre))
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = grouping_based_on_score) 

       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 PDSS.SR_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 PDSS.SR_2 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 PDSS.SR_3 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 PDSS.SR_4 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
5 PDSS.SR_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
6 PDSS.SR_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA

DIF-sex

Code

dif.sex <- dif %>% filter(Time=='PRE') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_pre, dif.sex)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.sex) 

       Item     Var   gamma     se pvalue padj.BH sig   lower  upper
1 PDSS.SR_1 dif.sex  0.0540 0.0534 0.3118  1.0000     -0.0507 0.1587
2 PDSS.SR_2 dif.sex  0.0513 0.0526 0.3290  1.0000     -0.0517 0.1543
3 PDSS.SR_3 dif.sex -0.0725 0.0496 0.1437  0.8625     -0.1698 0.0247
4 PDSS.SR_4 dif.sex -0.0911 0.0489 0.0627  0.3764     -0.1869 0.0048
5 PDSS.SR_6 dif.sex  0.0572 0.0484 0.2374  1.0000     -0.0377 0.1522
6 PDSS.SR_7 dif.sex  0.0195 0.0493 0.6924  1.0000     -0.0772 0.1162

DIF-age

Code

dif.age <- dif %>% filter(Time=='PRE') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_pre, dif.age)

Item	2	3	Mean location	StDev	MaxDiff
PDSS.SR_1	0.004	0.124	0.064	0.085	0.120
PDSS.SR_2	0.008	0.175	0.092	0.118	0.167
PDSS.SR_3	-1.030	-0.574	-0.802	0.322	0.456
PDSS.SR_4	-0.032	-0.362	-0.197	0.234	0.330
PDSS.SR_6	0.634	0.390	0.512	0.172	0.244
PDSS.SR_7	0.416	0.247	0.332	0.119	0.168

Code

partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.age) 

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 PDSS.SR_1 dif.age -0.0643 0.0316 0.0414  0.2486      -0.1262 -0.0025
2 PDSS.SR_2 dif.age -0.0418 0.0310 0.1780  1.0000      -0.1026  0.0190
3 PDSS.SR_3 dif.age -0.1211 0.0280 0.0000  0.0001  *** -0.1760 -0.0662
4 PDSS.SR_4 dif.age  0.1013 0.0292 0.0005  0.0031   **  0.0441  0.1584
5 PDSS.SR_6 dif.age  0.0354 0.0278 0.2023  1.0000      -0.0190  0.0899
6 PDSS.SR_7 dif.age  0.0485 0.0287 0.0907  0.5439      -0.0077  0.1047

DIF-education

Code

dif.edu <- dif %>% filter(Time=='PRE') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_pre, dif.edu)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.edu) 

       Item     Var   gamma     se pvalue padj.BH sig   lower   upper
1 PDSS.SR_1 dif.edu  0.1017 0.0473 0.0316  0.1894      0.0090  0.1945
2 PDSS.SR_2 dif.edu  0.0683 0.0471 0.1467  0.8800     -0.0239  0.1606
3 PDSS.SR_3 dif.edu  0.0585 0.0438 0.1815  1.0000     -0.0273  0.1442
4 PDSS.SR_4 dif.edu -0.0607 0.0436 0.1637  0.9821     -0.1462  0.0247
5 PDSS.SR_6 dif.edu -0.0150 0.0432 0.7288  1.0000     -0.0996  0.0696
6 PDSS.SR_7 dif.edu -0.1315 0.0438 0.0027  0.0162   * -0.2174 -0.0456

DIF-Year

Code

dif.yr <- dif %>% filter(Time=='PRE') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_pre, dif.yr)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.yr) 

       Item    Var   gamma     se pvalue padj.BH sig   lower  upper
1 PDSS.SR_1 dif.yr  0.0317 0.0330 0.3369  1.0000     -0.0330 0.0963
2 PDSS.SR_2 dif.yr -0.0090 0.0337 0.7901  1.0000     -0.0750 0.0571
3 PDSS.SR_3 dif.yr  0.0307 0.0301 0.3069  1.0000     -0.0282 0.0896
4 PDSS.SR_4 dif.yr -0.0525 0.0306 0.0864  0.5184     -0.1126 0.0075
5 PDSS.SR_6 dif.yr  0.0427 0.0302 0.1569  0.9414     -0.0164 0.1019
6 PDSS.SR_7 dif.yr -0.0357 0.0305 0.2411  1.0000     -0.0954 0.0240

9.1 Analysis part 2 decision

Item 2 displays the most misfit and most residual correlations (with 1 too).

10 Rasch analysis 3

Based on the above we remove item 2

Code

df_pre <- df %>% filter(Time=='PRE') %>% select(!c(Time,PDSS.SR_5,PDSS.SR_2)) #%>% sample_n(800)
RIlistItemsMargin(df_pre, fontsize = 12)

itemnr	item
PDSS.SR_1	Number
PDSS.SR_3	Worry/Fear
PDSS.SR_4	Avoidance places
PDSS.SR_6	Impairment work/home
PDSS.SR_7	Impairment social

Conditional item fit

Code

RIitemfit(df_pre, cutoff = "Smith98")

Item	InfitMSQ	Location	1 ± 2 / √n
PDSS.SR_1	1.151	0.04	[0.951, 1.049]
PDSS.SR_3	1.002	-0.84	[0.951, 1.049]
PDSS.SR_4	1.02	-0.11	[0.951, 1.049]
PDSS.SR_6	0.989	0.55	[0.951, 1.049]
PDSS.SR_7	0.869	0.36	[0.951, 1.049]
Note:
MSQ values based on conditional estimation (n = 1678 complete cases).

Code

simfit1 <- RIgetfit(df_pre, iterations = 300, cpu = n_cores)  

RIitemfit(df_pre, simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
PDSS.SR_1	1.151	[0.927, 1.088]	1.141	[0.928, 1.086]	0.063	0.055	0.04
PDSS.SR_3	1.002	[0.928, 1.069]	0.996	[0.931, 1.068]	no misfit	no misfit	-0.84
PDSS.SR_4	1.02	[0.935, 1.063]	1.023	[0.933, 1.075]	no misfit	no misfit	-0.11
PDSS.SR_6	0.989	[0.938, 1.077]	0.969	[0.931, 1.101]	no misfit	no misfit	0.55
PDSS.SR_7	0.869	[0.924, 1.054]	0.864	[0.922, 1.065]	0.055	0.058	0.36
Note:
MSQ values based on conditional calculations (n = 1678 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RIgetfitPlot(simfit1, df_pre)

CICC

Code

ICCplot(as.data.frame(df_pre), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[c(1,3,4,6)])

[1] "Please press Zoom on the Plots window to see the plot"

Code

### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)

Code

ICCplot(as.data.frame(df_pre), 
        itemnumber = 5, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[7])

Item-restscore

Code

RIrestscore(df_pre)

Item	Observed value	Model expected value	Absolute difference	Adjusted p-value (BH)	Statistical significance level	Location	Relative location
PDSS.SR_1	0.55	0.59	0.04	0.068	.	0.04	0.63
PDSS.SR_3	0.59	0.60	0.01	0.962		-0.84	-0.26
PDSS.SR_4	0.60	0.60	0.00	0.962		-0.11	0.47
PDSS.SR_6	0.61	0.60	0.01	0.883		0.55	1.14
PDSS.SR_7	0.67	0.60	0.07	0.000	***	0.36	0.94

Item-restscore-bootstrapped

Code

RIbootRestscore(df_pre,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600

Item	Item-restscore result	% of iterations	Conditional MSQ infit	Relative average item location
PDSS.SR_1	no misfit	88.4	1.15	0.63
PDSS.SR_1	underfit	11.6	1.15	0.63
PDSS.SR_3	no misfit	99.8	1.00	-0.26
PDSS.SR_4	no misfit	99.6	1.02	0.47
PDSS.SR_6	no misfit	97.8	0.99	1.14
PDSS.SR_7	no misfit	34.0	0.87	0.94
PDSS.SR_7	overfit	66.0	0.87	0.94
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 1678.

PCA

Code

RIpcmPCA(df_pre)

PCA of Rasch model residuals

Eigenvalues	Proportion of variance
1.55	32.2%
1.37	26.9%
1.06	21.7%
1.01	19%
0.01	0.1%

Residual correlations

Code

simcor1 <- RIgetResidCor(df_pre, iterations = 500, cpu = n_cores)
RIresidcorr(df_pre, cutoff = simcor1$p99)

	PDSS.SR_1	PDSS.SR_3	PDSS.SR_4	PDSS.SR_6	PDSS.SR_7
PDSS.SR_1
PDSS.SR_3	-0.07
PDSS.SR_4	-0.35	-0.17
PDSS.SR_6	-0.15	-0.28	-0.28
PDSS.SR_7	-0.29	-0.27	-0.03	-0.13
Note:
Relative cut-off value is -0.116, which is 0.085 above the average correlation (-0.201). Correlations above the cut-off are highlighted in red text.

Partial gamma LD

Code

RIpartgamLD(df_pre)

Item 1	Item 2	Partial gamma	SE	Lower CI	Upper CI	Adjusted p-value (BH)
PDSS.SR_7	PDSS.SR_4	0.314	0.039	0.237	0.391	0.000
PDSS.SR_4	PDSS.SR_7	0.224	0.041	0.144	0.304	0.000
PDSS.SR_3	PDSS.SR_1	0.203	0.041	0.123	0.284	0.000
PDSS.SR_7	PDSS.SR_6	0.197	0.041	0.117	0.277	0.000
PDSS.SR_1	PDSS.SR_3	0.16	0.041	0.079	0.241	0.002
PDSS.SR_6	PDSS.SR_1	0.138	0.042	0.055	0.220	0.021

1st contrast loadings

Code

RIloadLoc(df_pre)

Response categories

Code

mirt(df_pre, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code

# for fewer items or a more magnified figure, use:
#RIitemCats(df)

Targeting

Code

# increase fig-height above as needed, if you have many items
RItargeting(df_pre)

Item hierarchy

Code

RIitemHierarchy(df_pre)

Score groups

Code

iarm::score_groups(as.data.frame(df_pre)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)

  score_group   n
1           1 896
2           2 782

Code

dif_plots <- df_pre %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code

clr_tests(df_pre, model = "PCM")

Conditional Likelihood Ratio Tests:

        clr  df pvalue  sig 
overall 56.7 19 1.3e-05  ***

Code

item_obsexp(PCM(df_pre))

Score group 1: 
          mean obs mean exp std.res sig
PDSS.SR_1  1.170    1.146    1.283     
PDSS.SR_3  1.586    1.570    0.768     
PDSS.SR_4  1.207    1.212   -0.211     
PDSS.SR_6  0.761    0.752    0.478     
PDSS.SR_7  0.882    0.919   -1.829     

Score group 2: 
          mean obs mean exp std.res sig
PDSS.SR_1  2.082    2.109   -1.198     
PDSS.SR_3  2.704    2.714   -0.471     
PDSS.SR_4  2.388    2.383    0.215     
PDSS.SR_6  2.069    2.080   -0.419     
PDSS.SR_7  2.268    2.225    1.761     

Code

grouping_based_on_score = score_groups(as.data.frame(df_pre))
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = grouping_based_on_score) 

       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 PDSS.SR_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 PDSS.SR_3 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 PDSS.SR_4 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 PDSS.SR_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
5 PDSS.SR_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA

DIF-sex

Code

dif.sex <- dif %>% filter(Time=='PRE') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_pre, dif.sex)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.sex) 

       Item     Var   gamma     se pvalue padj.BH sig   lower  upper
1 PDSS.SR_1 dif.sex  0.0469 0.0506 0.3544  1.0000     -0.0523 0.1461
2 PDSS.SR_3 dif.sex -0.0569 0.0494 0.2490  1.0000     -0.1537 0.0398
3 PDSS.SR_4 dif.sex -0.0863 0.0496 0.0818  0.4089     -0.1835 0.0109
4 PDSS.SR_6 dif.sex  0.0618 0.0489 0.2057  1.0000     -0.0339 0.1576
5 PDSS.SR_7 dif.sex  0.0364 0.0505 0.4707  1.0000     -0.0625 0.1353

DIF-age

Code

dif.age <- dif %>% filter(Time=='PRE') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_pre, dif.age)

Item	2	3	Mean location	StDev	MaxDiff
PDSS.SR_1	0.000	0.149	0.075	0.105	0.149
PDSS.SR_3	-1.003	-0.525	-0.764	0.338	0.478
PDSS.SR_4	-0.028	-0.313	-0.171	0.201	0.285
PDSS.SR_6	0.621	0.410	0.516	0.149	0.211
PDSS.SR_7	0.409	0.278	0.344	0.093	0.131

Code

partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.age) 

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 PDSS.SR_1 dif.age -0.0710 0.0302 0.0187  0.0933    . -0.1301 -0.0118
2 PDSS.SR_3 dif.age -0.1377 0.0277 0.0000  0.0000  *** -0.1920 -0.0834
3 PDSS.SR_4 dif.age  0.0991 0.0294 0.0008  0.0038   **  0.0414  0.1568
4 PDSS.SR_6 dif.age  0.0393 0.0284 0.1666  0.8332      -0.0164  0.0950
5 PDSS.SR_7 dif.age  0.0446 0.0294 0.1297  0.6483      -0.0131  0.1023

DIF-education

Code

dif.edu <- dif %>% filter(Time=='PRE') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_pre, dif.edu)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.edu) 

       Item     Var   gamma     se pvalue padj.BH sig   lower   upper
1 PDSS.SR_1 dif.edu  0.0824 0.0457 0.0711  0.3557     -0.0071  0.1719
2 PDSS.SR_3 dif.edu  0.0803 0.0438 0.0669  0.3344     -0.0056  0.1661
3 PDSS.SR_4 dif.edu -0.0374 0.0439 0.3942  1.0000     -0.1235  0.0486
4 PDSS.SR_6 dif.edu -0.0055 0.0441 0.9007  1.0000     -0.0920  0.0810
5 PDSS.SR_7 dif.edu -0.1295 0.0450 0.0040  0.0198   * -0.2176 -0.0414

DIF-Year

Code

dif.yr <- dif %>% filter(Time=='PRE') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_pre, dif.yr)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.yr) 

       Item    Var   gamma     se pvalue padj.BH sig   lower  upper
1 PDSS.SR_1 dif.yr  0.0310 0.0319 0.3319  1.0000     -0.0316 0.0935
2 PDSS.SR_3 dif.yr  0.0422 0.0297 0.1551  0.7754     -0.0160 0.1004
3 PDSS.SR_4 dif.yr -0.0599 0.0313 0.0552  0.2760     -0.1212 0.0013
4 PDSS.SR_6 dif.yr  0.0414 0.0305 0.1739  0.8695     -0.0183 0.1011
5 PDSS.SR_7 dif.yr -0.0573 0.0308 0.0626  0.3129     -0.1176 0.0030

10.1 Analysis part 3 decision

The only really underperforming item is item 7. With the overfit indicated by the boostrap and the association / local dependence between item 4 and 7 as indicated by the partial gamma.

11 Rasch analysis 4

Based on the above we remove item 7

Code

df_pre <- df %>% filter(Time=='PRE') %>% select(!c(Time,PDSS.SR_5,PDSS.SR_2,PDSS.SR_7)) #%>% sample_n(800)
RIlistItemsMargin(df_pre, fontsize = 12)

itemnr	item
PDSS.SR_1	Number
PDSS.SR_3	Worry/Fear
PDSS.SR_4	Avoidance places
PDSS.SR_6	Impairment work/home

Conditional item fit

Code

RIitemfit(df_pre, cutoff = "Smith98")

Item	InfitMSQ	Location	1 ± 2 / √n
PDSS.SR_1	1.053	0.12	[0.951, 1.049]
PDSS.SR_3	0.92	-0.72	[0.951, 1.049]
PDSS.SR_4	1.07	-0.02	[0.951, 1.049]
PDSS.SR_6	0.993	0.62	[0.951, 1.049]
Note:
MSQ values based on conditional estimation (n = 1678 complete cases).

Code

simfit1 <- RIgetfit(df_pre, iterations = 300, cpu = n_cores)  

RIitemfit(df_pre, simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
PDSS.SR_1	1.053	[0.916, 1.076]	1.043	[0.921, 1.077]	no misfit	no misfit	0.12
PDSS.SR_3	0.92	[0.931, 1.072]	0.918	[0.93, 1.077]	0.011	0.012	-0.72
PDSS.SR_4	1.07	[0.929, 1.074]	1.077	[0.927, 1.08]	no misfit	no misfit	-0.02
PDSS.SR_6	0.993	[0.934, 1.064]	0.974	[0.915, 1.084]	no misfit	no misfit	0.62
Note:
MSQ values based on conditional calculations (n = 1678 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RIgetfitPlot(simfit1, df_pre)

CICC

Code

ICCplot(as.data.frame(df_pre), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[c(1,3,4,6)])

[1] "Please press Zoom on the Plots window to see the plot"

Code

### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)

Item-restscore

Code

RIrestscore(df_pre)

Item	Observed value	Model expected value	Absolute difference	Adjusted p-value (BH)	Statistical significance level	Location	Relative location
PDSS.SR_1	0.56	0.56	0.00	0.748		0.12	0.59
PDSS.SR_3	0.60	0.56	0.04	0.049	*	-0.72	-0.25
PDSS.SR_4	0.54	0.56	0.02	0.462		-0.02	0.45
PDSS.SR_6	0.57	0.57	0.00	0.748		0.62	1.09

Item-restscore-bootstrapped

Code

RIbootRestscore(df_pre,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600

Item	Item-restscore result	% of iterations	Conditional MSQ infit	Relative average item location
PDSS.SR_1	no misfit	99.8	1.05	0.59
PDSS.SR_3	no misfit	83.8	0.92	-0.25
PDSS.SR_3	overfit	16.2	0.92	-0.25
PDSS.SR_4	no misfit	99.8	1.07	0.45
PDSS.SR_6	no misfit	100.0	0.99	1.09
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 1678.

PCA

Code

RIpcmPCA(df_pre)

PCA of Rasch model residuals

Eigenvalues	Proportion of variance
1.47	38.8%
1.41	33.4%
1.11	27.6%
0.01	0.2%

Residual correlations

Code

simcor1 <- RIgetResidCor(df_pre, iterations = 500, cpu = n_cores)
RIresidcorr(df_pre, cutoff = simcor1$p99)

	PDSS.SR_1	PDSS.SR_3	PDSS.SR_4	PDSS.SR_6
PDSS.SR_1
PDSS.SR_3	-0.14
PDSS.SR_4	-0.37	-0.18
PDSS.SR_6	-0.19	-0.32	-0.25
Note:
Relative cut-off value is -0.162, which is 0.081 above the average correlation (-0.242). Correlations above the cut-off are highlighted in red text.

Partial gamma LD

Code

RIpartgamLD(df_pre)

Item 1	Item 2	Partial gamma	SE	Lower CI	Upper CI	Adjusted p-value (BH)
PDSS.SR_3	PDSS.SR_4	0.175	0.042	0.093	0.257	0.000
PDSS.SR_3	PDSS.SR_1	0.152	0.044	0.065	0.238	0.007
PDSS.SR_6	PDSS.SR_1	0.13	0.043	0.045	0.215	0.031
PDSS.SR_1	PDSS.SR_6	0.122	0.042	0.039	0.205	0.049

1st contrast loadings

Code

RIloadLoc(df_pre)

Response categories

Code

mirt(df_pre, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code

# for fewer items or a more magnified figure, use:
#RIitemCats(df)

Targeting

Code

# increase fig-height above as needed, if you have many items
RItargeting(df_pre)

Item hierarchy

Code

RIitemHierarchy(df_pre)

Score groups

Code

iarm::score_groups(as.data.frame(df_pre)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)

  score_group   n
1           1 790
2           2 888

Code

dif_plots <- df_pre %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code

clr_tests(df_pre, model = "PCM")

Conditional Likelihood Ratio Tests:

        clr  df pvalue  sig 
overall 38.5 15 0.00076  ***

Code

item_obsexp(PCM(df_pre))

Score group 1: 
          mean obs mean exp std.res sig
PDSS.SR_1  1.069    1.089   -1.002     
PDSS.SR_3  1.485    1.493   -0.340     
PDSS.SR_4  1.150    1.131    0.864     
PDSS.SR_6  0.684    0.672    0.564     

Score group 2: 
          mean obs mean exp std.res sig
PDSS.SR_1  2.061    2.044    0.832     
PDSS.SR_3  2.659    2.650    0.441     
PDSS.SR_4  2.297    2.314   -0.764     
PDSS.SR_6  1.980    1.990   -0.437     

Code

grouping_based_on_score = score_groups(as.data.frame(df_pre))
partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = grouping_based_on_score) 

       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 PDSS.SR_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 PDSS.SR_3 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 PDSS.SR_4 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 PDSS.SR_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA

DIF-sex

Code

dif.sex <- dif %>% filter(Time=='PRE') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_pre, dif.sex)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.sex) 

       Item     Var   gamma     se pvalue padj.BH sig   lower  upper
1 PDSS.SR_1 dif.sex  0.0694 0.0522 0.1841  0.7364     -0.0330 0.1718
2 PDSS.SR_3 dif.sex -0.0514 0.0511 0.3141  1.0000     -0.1516 0.0487
3 PDSS.SR_4 dif.sex -0.0816 0.0492 0.0969  0.3877     -0.1780 0.0147
4 PDSS.SR_6 dif.sex  0.0723 0.0489 0.1398  0.5592     -0.0237 0.1682

DIF-age

Code

dif.age <- dif %>% filter(Time=='PRE') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_pre, dif.age)

Item	2	3	Mean location	StDev	MaxDiff
PDSS.SR_1	0.093	0.201	0.147	0.077	0.109
PDSS.SR_3	-0.867	-0.430	-0.648	0.309	0.437
PDSS.SR_4	0.073	-0.227	-0.077	0.212	0.299
PDSS.SR_6	0.701	0.455	0.578	0.174	0.246

Code

partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.age) 

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 PDSS.SR_1 dif.age -0.0643 0.0312 0.0391  0.1563      -0.1253 -0.0032
2 PDSS.SR_3 dif.age -0.1356 0.0291 0.0000  0.0000  *** -0.1927 -0.0786
3 PDSS.SR_4 dif.age  0.1023 0.0287 0.0004  0.0014   **  0.0462  0.1585
4 PDSS.SR_6 dif.age  0.0487 0.0286 0.0885  0.3539      -0.0073  0.1047

DIF-education

Code

dif.edu <- dif %>% filter(Time=='PRE') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_pre, dif.edu)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.edu) 

       Item     Var   gamma     se pvalue padj.BH sig   lower  upper
1 PDSS.SR_1 dif.edu  0.0600 0.0469 0.2011  0.8042     -0.0320 0.1519
2 PDSS.SR_3 dif.edu  0.0518 0.0452 0.2511  1.0000     -0.0367 0.1403
3 PDSS.SR_4 dif.edu -0.0803 0.0434 0.0646  0.2584     -0.1654 0.0049
4 PDSS.SR_6 dif.edu -0.0279 0.0440 0.5259  1.0000     -0.1141 0.0583

DIF-Year

Code

dif.yr <- dif %>% filter(Time=='PRE') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_pre, dif.yr)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_pre),
            dat.exo = dif.yr) 

       Item    Var   gamma     se pvalue padj.BH sig   lower   upper
1 PDSS.SR_1 dif.yr  0.0158 0.0325 0.6261  1.0000     -0.0478  0.0794
2 PDSS.SR_3 dif.yr  0.0201 0.0311 0.5185  1.0000     -0.0409  0.0811
3 PDSS.SR_4 dif.yr -0.0661 0.0308 0.0321  0.1285     -0.1265 -0.0056
4 PDSS.SR_6 dif.yr  0.0277 0.0309 0.3688  1.0000     -0.0328  0.0883

11.1 Analysis part 4 decision

While item 3 have some small overfit in the boostrap (~15%) it is not so severe. As is the small DIF indications on age.

12 Analysis conclusion

Resulting PDSS-SR reworked scale thus has item 1, 3, 4 and 6.

12.1 Resulting items

Code

df_pre <- df %>% filter(Time=='PRE') %>% select(!c(Time,PDSS.SR_5,PDSS.SR_2,PDSS.SR_7)) #%>% sample_n(800)

what_scale <- gsub('\\.','',items_to_use)

Code

RIlistItemsMargin(df_pre, fontsize = 12)

itemnr	item
PDSS.SR_1	Number
PDSS.SR_3	Worry/Fear
PDSS.SR_4	Avoidance places
PDSS.SR_6	Impairment work/home

12.2 Reliability

Code

RItif(df_pre)

Code

RItif(df_pre,samplePSI=TRUE)

12.3 Person parameter

Code

RIpfit(df_pre)

12.4 Misfit check

Code

pfit_u3poly <- PerFit::U3poly(matrix = df_pre, 
                      Ncat = 5, # make sure to input number of response categories, not thresholds
                      IRT.PModel = "PCM")
misfits <- PerFit::flagged.resp(pfit_u3poly) %>% 
  pluck("Scores") %>% 
  as.data.frame() %>% 
  pull(FlaggedID)

misfits2 <- RIpfit(df_pre, output = "rowid")

All resp

Code

RIitemfit(df_pre, simcut = simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
PDSS.SR_1	1.053	[0.916, 1.076]	1.043	[0.921, 1.077]	no misfit	no misfit	0.12
PDSS.SR_3	0.92	[0.931, 1.072]	0.918	[0.93, 1.077]	0.011	0.012	-0.72
PDSS.SR_4	1.07	[0.929, 1.074]	1.077	[0.927, 1.08]	no misfit	no misfit	-0.02
PDSS.SR_6	0.993	[0.934, 1.064]	0.974	[0.915, 1.084]	no misfit	no misfit	0.62
Note:
MSQ values based on conditional calculations (n = 1678 complete cases). Simulation based thresholds from 300 simulated datasets.

U3 misfit removed

Code

RIitemfit(df_pre[-misfits,], simcut = simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
PDSS.SR_1	1.006	[0.916, 1.076]	0.991	[0.921, 1.077]	no misfit	no misfit	0.97
PDSS.SR_3	0.98	[0.931, 1.072]	0.937	[0.93, 1.077]	no misfit	no misfit	-3.50
PDSS.SR_4	1.017	[0.929, 1.074]	1.012	[0.927, 1.08]	no misfit	no misfit	0.79
PDSS.SR_6	1.015	[0.934, 1.064]	0.977	[0.915, 1.084]	no misfit	no misfit	1.74
Note:
MSQ values based on conditional calculations (n = 1495 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RItif(df_pre[-misfits,])

ZSTD removed

Code

RIitemfit(df_pre[-misfits2,], simcut = simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
PDSS.SR_1	1.073	[0.916, 1.076]	1.065	[0.921, 1.077]	no misfit	no misfit	0.15
PDSS.SR_3	0.942	[0.931, 1.072]	0.935	[0.93, 1.077]	no misfit	no misfit	-0.79
PDSS.SR_4	1.016	[0.929, 1.074]	1.023	[0.927, 1.08]	no misfit	no misfit	-0.02
PDSS.SR_6	0.998	[0.934, 1.064]	0.978	[0.915, 1.084]	no misfit	no misfit	0.66
Note:
MSQ values based on conditional calculations (n = 1559 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RItif(df_pre[-misfits2,])

12.5 Item parameters

Code

RIitemparams(df_pre)

	Threshold 1	Threshold 2	Threshold 3	Threshold 4	Item location
PDSS.SR_1	-3.17	-0.54	1.88	2.30	0.12
PDSS.SR_3	-4.02	-1.40	0.11	2.42	-0.72
PDSS.SR_4	-2.70	-0.93	0.97	2.59	-0.02
PDSS.SR_6	-1.69	0.04	1.05	3.07	0.62
Note:
Item location is the average of the thresholds for each item.

Code

item_param_matrix <- RIitemparams(df_pre,output='dataframe') %>% as.matrix(.)
write.csv(item_param_matrix,file=paste0("./results/item_params_",what_scale,'.csv'),row.names=TRUE)

12.6 Ordinal sum to interval score

Code

RIscoreSE(df_pre)

Ordinal sum score	Logit score	Logit std.error
0	-5.419	0.775
1	-3.960	0.968
2	-3.072	0.901
3	-2.367	0.818
4	-1.770	0.760
5	-1.243	0.721
6	-0.759	0.696
7	-0.300	0.679
8	0.146	0.669
9	0.587	0.663
10	1.025	0.662
11	1.460	0.668
12	1.892	0.689
13	2.340	0.731
14	2.849	0.793
15	3.523	0.842
16	4.787	0.693

Code

RIscoreSE(df_pre,output='figure')

Code

sum_to_latent <- RIscoreSE(df_pre,output = 'dataframe')
write.csv(sum_to_latent,file=paste0('./results/ordinal_sum_to_latent_',what_scale,'.csv'),row.names=FALSE)

12.7 Thetas

The ordinal sum to interval score contains the information to transform into the below thetas, but we plot them here for convinience (based on 4 items)

Code

est_thetas2 <- RIestThetas(df_pre, method = "WLE")
hist(est_thetas2$WLE, 
     col = "#009ca6", 
     main = "Histogram of person locations (thetas)", 
     breaks = 10)

13 Software used

Code

pkgs <- grateful::cite_packages(cite.tidyverse = TRUE, 
                      output = "table",
                      bib.file = "grateful-refs.bib",
                      include.RStudio = FALSE,
                      out.dir = getwd())
formattable(pkgs, table.attr = 'class=\"table table-striped\" style="font-size: 15px; font-family: Lato; width: 80%"')

Package	Version	Citation
base	4.4.1	@base
car	3.1.3	@car
doParallel	1.0.17	@doParallel
easyRasch	0.3.3	@easyRasch
eRm	1.0.6	@eRm2007a; @eRm2007b; @eRm2009c; @eRm2013d; @eRm2015e; @eRm2019f
formattable	0.2.1	@formattable
furrr	0.3.1	@furrr
ggrepel	0.9.6	@ggrepel
glue	1.8.0	@glue
gridExtra	2.3	@gridExtra
iarm	0.4.3	@iarm
janitor	2.2.0	@janitor
kableExtra	1.4.0	@kableExtra
knitr	1.49	@knitr2014; @knitr2015; @knitr2024
matrixStats	1.4.1	@matrixStats
mirt	1.43	@mirt
mokken	3.1.2	@mokken2007; @mokken2012
patchwork	1.3.0	@patchwork
PerFit	1.4.6	@PerFit
psych	2.4.6.26	@psych
psychotools	0.7.4	@psychotools2021; @psychotools2022; @psychotools2024
psychotree	0.16.1	@psychotree2010e; @psychotree2011a; @psychotree2015b; @psychotree2018c; @psychotree2018d
reshape	0.8.9	@reshape
rmarkdown	2.29	@rmarkdown2018; @rmarkdown2020; @rmarkdown2024
tidyverse	2.0.0	@tidyverse

14 Session info

Code

sessionInfo()

R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Stockholm
tzcode source: system (glibc)

attached base packages:
 [1] parallel  grid      stats4    stats     graphics  grDevices utils    
 [8] datasets  methods   base     

other attached packages:
 [1] gridExtra_2.3     knitr_1.49        car_3.1-3         carData_3.0-5    
 [5] grateful_0.2.10   easyRasch_0.3.3   doParallel_1.0.17 iterators_1.0.14 
 [9] furrr_0.3.1       future_1.34.0     foreach_1.5.2     ggdist_3.3.2     
[13] janitor_2.2.0     iarm_0.4.3        hexbin_1.28.5     catR_3.17        
[17] glue_1.8.0        ggrepel_0.9.6     patchwork_1.3.0   reshape_0.8.9    
[21] matrixStats_1.4.1 psychotree_0.16-1 psychotools_0.7-4 partykit_1.2-23  
[25] mvtnorm_1.3-2     libcoin_1.0-10    psych_2.4.6.26    mirt_1.43        
[29] lattice_0.22-5    eRm_1.0-6         lubridate_1.9.4   forcats_1.0.0    
[33] stringr_1.5.1     dplyr_1.1.4       purrr_1.0.2       readr_2.1.5      
[37] tidyr_1.3.1       tibble_3.2.1      ggplot2_3.5.1     tidyverse_2.0.0  
[41] kableExtra_1.4.0  formattable_0.2.1

loaded via a namespace (and not attached):
  [1] splines_4.4.1        bitops_1.0-6         R.oo_1.27.0         
  [4] cellranger_1.1.0     deSolve_1.40         rpart_4.1.23        
  [7] lifecycle_1.0.4      rprojroot_2.0.4      globals_0.16.3      
 [10] MASS_7.3-61          backports_1.5.0      magrittr_2.0.3      
 [13] vcd_1.4-13           Hmisc_5.2-1          rmarkdown_2.29      
 [16] yaml_2.3.10          sm_2.2-6.0           sessioninfo_1.2.2   
 [19] pbapply_1.7-2        RColorBrewer_1.1-3   abind_1.4-8         
 [22] audio_0.1-11         expm_1.0-0           quadprog_1.5-8      
 [25] R.utils_2.12.3       RCurl_1.98-1.16      pracma_2.4.4        
 [28] nnet_7.3-19          listenv_0.9.1        testthat_3.2.2      
 [31] RPushbullet_0.3.4    vegan_2.6-8          parallelly_1.40.1   
 [34] svglite_2.1.3        permute_0.9-7        codetools_0.2-19    
 [37] xml2_1.3.6           tidyselect_1.2.1     farver_2.1.2        
 [40] base64enc_0.1-3      jsonlite_1.8.9       polycor_0.8-1       
 [43] ks_1.14.3            progressr_0.15.1     Formula_1.2-5       
 [46] survival_3.7-0       systemfonts_1.1.0    tools_4.4.1         
 [49] gnm_1.1-5            ragg_1.3.3           snow_0.4-4          
 [52] Rcpp_1.0.13-1        mnormt_2.1.1         admisc_0.36         
 [55] xfun_0.49            here_1.0.1           mgcv_1.9-1          
 [58] msm_1.8.2            distributional_0.5.0 ca_0.71.1           
 [61] withr_3.0.2          beepr_2.0            fastmap_1.2.0       
 [64] fansi_1.0.6          digest_0.6.37        timechange_0.3.0    
 [67] R6_2.5.1             textshaping_0.4.1    colorspace_2.1-1    
 [70] R.methodsS3_1.8.2    inum_1.0-5           utf8_1.2.4          
 [73] generics_0.1.3       renv_1.0.10          data.table_1.16.4   
 [76] SimDesign_2.17.1     htmlwidgets_1.6.4    pkgconfig_2.0.3     
 [79] gtable_0.3.6         lmtest_0.9-40        pcaPP_2.0-5         
 [82] brio_1.1.5           htmltools_0.5.8.1    scales_1.3.0        
 [85] snakecase_0.11.1     irtoys_0.2.2         rstudioapi_0.17.1   
 [88] tzdb_0.4.0           checkmate_2.3.2      nlme_3.1-165        
 [91] curl_6.0.1           zoo_1.8-12           KernSmooth_2.23-24  
 [94] relimp_1.0-5         vcdExtra_0.8-5       fda_6.2.0           
 [97] foreign_0.8-86       pillar_1.9.0         vctrs_0.6.5         
[100] Deriv_4.1.6          cluster_2.1.6        dcurver_0.9.2       
[103] GPArotation_2024.3-1 htmlTable_2.4.3      evaluate_1.0.1      
[106] cli_3.6.3            compiler_4.4.1       rlang_1.1.4         
[109] future.apply_1.11.3  labeling_0.4.3       mclust_6.1.1        
[112] fds_1.8              plyr_1.8.9           rainbow_3.8         
[115] stringi_1.8.4        hdrcde_3.4           viridisLite_0.4.2   
[118] munsell_0.5.1        Matrix_1.6-5         qvcalc_1.0.3        
[121] hms_1.1.3            ltm_1.2-0            PerFit_1.4.6        
[124] readxl_1.4.3        

15 References

Reuse

https://creativecommons.org/licenses/by/4.0/